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Abstract 

It is known that Goertzel's algorithm is much less numerically 
accurate than the Fast Fourier Transform (FFT)(Cf. 2 ). In order 
to improve accuracy we propose modifications of both Goertzel's and 
Horner's algorithms based on the divide-and-conquer techniques. The 
proof of the numerical stability of these two modified algorithms is 
given. The numerical tests in Matlab demonstrate the computational 
advantages of the proposed modifications. The appendix contains 
the proof of numerical stability of Goertzel's algorithm of polynomial 
evaluation. 

AMS subject classification: 65F35, 65G50. 
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1 Introduction 

The aim of this paper is to improve the accm'acy of polynomial evaluation, 
mainly Horner's and Goertzel's algorithms. Both, Horner's and Goertzel's 
methods are frequently used in the interpolation and approximation prob- 
lems and in signal processing. Goertzel's algorithm is implemented in Mat- 
lab, it's included in the Signal Processing Toolbox. The function "fft" re- 
turns the Discrete Fourier Transform (DFT) computed with a Fast Fourier 
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Transform (FFT) algorithm and the function "goertzel" computes DFT of 
specific indices in a vector. 
In this paper we consider more general case of evaluating a polynomial 

N 

(1) w{z) = ^anz'', 

n=0 

where z £ C and ao, • • • , cat S C. 
Note that for z = e*^ and oq, • • • , cin-,^ € M we have 

N N 

w{e^^) = an cos nS, + i a„ sin n^. 

n=0 n=0 

DFT returns yk = w{zk), k = 0,. . . ,N, where Zk = S'' = (Cf. 0, 

p. 10). 

It is observed (see "help goertzel" in Matlab Signal Processing Tool- 
box) that compared with the Fast Fourier Transform algorithm (FFT), 
Goertzel's algorithm is much less numerically accurate, which can be visi- 
ble especially for high-scale problems. 

We propose the algorithm PEMA (Polynomial Evaluation Modified Al- 
gorithm), which is based on the repetitive use of some algorithm W for 
evaluating polynomials. This algorithm can be e.g. Horner's or Goertzel's 
scheme. The cost of PEMA is comparable to the cost of W and the error 
bound of PEMA may be significantly smaller than the error bound of W . 
We prove that if W is stable then PEMA is also numerically stable (see 
section 3.2). In practice, one should use only numerically stable algorithms. 

We say that an algorithm of evaluating (1) is componentwise backward 
stable with respect to the data oq, • • • , aat G C and z G C if the value w{z) 
computed by this algorithm is an exact value of a polynomial for slightly 
perturbed coefficients a„ and z, i.e. 

N 

(2) w{z) = ^[an{l + l^n)\[z{l + (3)T, \^^n\<ANeM, \P\<ZNeM, 

n=0 

where and Zj\f are modestly growing functions of and eM is the 
machine precision. 



Horner and Goertzel 



3 



Throughout the paper we assume that the coefficients of a polynomial 
w{z) are complex. 

In the error analysis of PEMA we consider perturbations not only of 
polynomial coefficients, but also of z. Notice that usually the exact value 
of z is not known, e.g. z is given as z = e*^. Then z = c + is, c = cos^, 
s = sin^ and c = c + Ac, s = s + As, |Ac|, |As| < I'eM, where u is small. 
Then the perturbed value z can be written as z = z (1 + ry), \ri\ < \plv eu- 

Then with help of Taylor expansion (2) leads to 

AT 

\w{z)-w{z)\<\(i\\zw'{z)\+Y.\an\\bn\\z\' ^0{el,), 

n=0 

and further 

N 

\w{z)-w{z)\< eA/(^ivX] +^^1-^1 l^'WI) + C>i4i)- 

n=0 

Numerical stability of Horner's algorithm was first given by Wilkinson 
(Cf. [13], pp. 36-37, 49-50) who proved that Zn = and An ~ 2N, pro- 
vided that the data gq, . . . , S K and z G M are exactly representable in 
floating point arithmetic (fi). Despite of a bad reputation of Goertzel's 
algorithm as a method of computing Fourier series '}2in=o cos and 
Yln=i sinn^ with respect to the data oq, • • • , oat G M and a given ar- 
gument C G IK (Cf. [m, pp. 84-88, H, m. El) we prove that Goertzel's 
algorithm is numerically stable in a sense (2). The respective constants are 
Zn = and An is of order A^^, provided that the data oq, . . ■ , oa? and z 
are exactly representable in fl (see Theorem 2 and Table 0). 

In order to improve accuracy we propose modifications of both Goertzel's 
and Horner's algorithms based on the divide-and-conquer techniques. The 
idea is not quite new, there are numerous divide-and-conquer parallel al- 
gorithms for polynomial evaluation (Cf. [Hj, p. 70). The goal of our work 
is to split a polynomial in "the proper way" in order to refine results. We 
show that the constants An and Zn in (2) can be significantly decreased, 
in comparison with the classical Horner's and Goertzel's algorithms, which 
is of great importance for large N (see Table in section 3), e.g. for N = 2^ 
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our divide-and-conquer algorithm PEMA results in A]sf of order log2N and 
Zn of order unity. 

Tests included in section 4 confirm theoretical results. We also imple- 
mented Reinsch's modification of Goertzel's algorithm (Cf. 11 , pp. 86-88) 
for evaluation of (1), but it turned out that the numerical results were com- 
parable to these given by standard Goertzel's algorithm. For this reason 
we don't include them in section 4 devoted to numerical experiments. 

2 Classical polynomial evaluation schemes 

The Horner scheme is the standard method for evaluation of a polynomial 
(1) at a given point z G C. We assume that oq, . . . , aiv G C. We write w{z) 
as follows 

w{z) = ao + z[ai + z(. . . + z{aN-i + zqn) ■ ■ ■))■ 

Algorithm 1 (Horner's rule) 

w:=0 

for n = N,N -1,...,0 

w := Un + zw 

end 

w{z) := w 

The complexity of Horner's algorithm Cn{H), counted as a number of 
multiplications is equal to N, which gives in general Cn{H) = 4:N real 
multiplications. We assume that the product of two complex numbers is 
computed in a natural way and in consequence one complex multiplication 
is equivalent to four real ones. 

The idea of Goertzel's algorithm is different. Suppose z = x + iy. Divide 

N 

a polynomial 'w{X) = ^ flnA" by a quadratic polynomial (A — z){X — z) = 

n=0 

X — pX — q with real coefficients p and q, where p = 2x and q = —\z\ . 
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Then 

N 

w{X) = {X-z){X-z)J2 ^nA"-^ + bo + biX 

n=2 

and, consequently, w{z) = 60 + 61 z. This leads to the following 
Algorithm 2 (Goertzel's algorithm) 

p := 2x 

q := -{x^+y^) 
bN+i •■= 
bN '■= Qn 
for n = N 

bn --an+p bn+l + q bn+2 

end 

u := (ao + xbi + qb2) 

v:=ybi 

w{z) := u + iv 

In general, the number of real multiplications needed by Goertzel's method 
is the same as those needed by Horner's algorithm. However, in special 
cases each of these algorithms can be less expensive than the other. For 
example, for z G M Goertzel's algorithm is twice as expensive as Horner's 
rule regardless of the polynomial coefficients. On the other hand consider 
the case of polynomial with real coefficients and 2: G C , = 1. Then all 
bn are real, q = 1 and bn = an + p 6n+i + for n = AT — 1, . . . , 1. 
The complexity of Goertzel's method reduces to N while the cost of Horner's 
rule is still 4N. 

Note that if z = 1, then w{z) = '^n=o'^n and Horner's rule is nothing 
else but a backward summation. 
We now derive an algorithm based on the divide-and-conquer technique. 



Horner and Goertzel 



6 



3 A new polynomial evaluation modified algorithm 
(PEMA) 

N 

Suppose a polynomial w{z) is given by iv{z) = ^2 o-nZ^ where N = and 

n=0 

s > 1. We can write ■w{z) in the following form: 

w{z) = {ao+aiz-\ \-as-iz^~^}+{as+as+iz-\ \-a2s-iz^~^}z^ -\ h 

+ {a(^sP-'^-i)s + a(sP-i-i)s+iZ-\ hasP-iz^^^jiz'^y" + asp{z''Y'' ' = 

„p-i 



j=0 

where zi = z^ af = Efe=o J = 0, 1> • • • , - 1, = 

and a^'^^ = aj for j = 0, 1, . . . , A/". 

Now we can interpret Z^jLo '^j^^^i ^ ^ polynomial of variable zi with the 
coefficients a^-^^ and proceed in the same manner as before. We continue 
this process and for m = 0, 1, . . . ,p — 1 write w{z) as follows 

gP-m 
3=0 

where zq = z and Zm = -^^-i m = 1,2, . . . ,p — 1. 
It is easy to prove that for m = 1, 2, . . . ,p — 1 and j = 0,1, . . . , s^""*— 1 



(3) Oj ^ — ^ ] ajsm_^_j. z^ . 

r=0 



For complexity and computational accuracy reasons we don't evaluate (3) 
directly, by Horner or Goertzel algorithm for polynomial of variable z and 
degree s'^ — 1, but use the relation 



s-l 

(m) ("^-1) k 



'j ~ ^js+k 
k=0 

Notice that Uj^^ is a polynomial of variable Zm-i and degree s — 1. 
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More precisely, given an algorithm W for evaluating polynomials, e.g. 
Horner's or Goertzel's algorithm, we produce a new divide-and-conquer 
algorithm. 

Algorithm 3 (PEMA) 

This algorithm uses the divide-and-conquer method to compute w{z) 
where z G M or z G C. The coefficients a„ may be either complex or 
real. 

1. zo = z 

af^ =aj for j = 0,l,...,iV 

2. for m = 1, . . . ,p — 1 

(m) (m— 1) 

flgp-m — agp_(„_i) — CLN 

for j = 0,l,...,sP-'"-l 

compute a^™'* = ^ o!f^,^\^_i by algorithm 
fc=o 

end 
end 

3. compute w{z) = ^ ^l-i algorithm W 

Note that p = 1 implies N = s and PEMA is nothing else but W applied 
to w^z). 

PEMA is an extension of a summation algorithm proposed in j4j. For 
N = 2^ and z = 1 PEMA coincides with the log-sum algorithm. 

3.1 Total cost of PEMA 

Suppose the complexity of the algorithm W is Cn = bN, b = const, i.e. W 

N 

needs Cat multiplications to compute w^z) = ^ dnZ^- We give a formula 

n=0 
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for complexity of PEMA valid under assumption that Zm is computed in a 
natural way, (see section 3.2): 



CiPEMA)=Y^{ ^ C7,_i+(s-l)}+C, = (s-l)(p-l)+C,+a_i^^^-^ 

m=l j=0 

According to this formula the complexity of PEMA with Horner is equal 
to Cat + (s — l){p — 1). Very often the latter term is not significant in 
comparison with Cat. 

Remark. Each a^j^\ j = 0, 1, . . . , s^^™ — 1 can be computed indepen- 
dently. It's a big advantage of PEMA because of possibility of parallel 
implementation, which can be useful especially for really large problems. 

3.2 Error analysis of PEMA 

We consider complex arithmetic (cfl) implemented using standard real 
arithmetic with machine precision e^,/- Then 

(4) cfl{x + y) = {x + y){l + 6), \ 6 \< cm foTx,yeC 

and provided that the product xy is computed using an ordinary algorithm 
we have (Cf. [5]) 

(5) cfl{xy) = {xy){l + r]), | r/ |< cem, 
where 



(6) 



1 for x,y eR or x G M, y € C 
1 + \/2 for x,yeC. 



The value Zm = z^-i is determined in a natural way by computing the 

consecutive powers of Zm, i.e. Zm-i, -Zm-i) ■ • • ; ^m-i- 

Then 

(7) z^ = c//(z^_i) = z^_i(l + 5^), \5m\<{s-l)ceM + 0{elj). 



Now we are in a position to give the error analysis of the PEMA algorithm. 
For simplicity we assume that oq, . . . , oat and z are represented exactly in 
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cfi and that s and p are fixed, N = s'^. We also assume that the result 
given by the algorithm W of evaluating w{z) = "^^=0 ^nz'^ in cfl satisfies 

N 

(8) =^ a„(l + A„)z", \An\<ANeM, 

n=0 

where An is an increasing function of N. W in PEMA can be Horner's or 
Goertzel's rule. For detailed information on Aj^ see (37). 

For m = I, . . . ,p—l and j = 0,1, . . . , s^~"* — 1 the values aj™"* , computed 
in cfl, can be written as follows 

(9) = X: a^^il + aJ.J ) 4-1, I aJ.J |< cm- 

k=0 

The formula (7) allows us to write Zm-i in the following way 

m— 1 

(10) Zm-1 = [Z{1 + ^m)r"'~\ l+Jm=ll{l + 6t)^. 

t=l 

Prom (7) we obtain an upper bound for | 7^ | 

m— 1 ^ 

I 7m |< (S - 1) CCM XI - + ^(^m)- 

t=l * 

Thus 

(11) \lm\<ceM + 0{ei^). 

Lemma 1 Assume that cfl{z) = z and cfl{an) = an, n = 0, . . . ,N and 
N = s'^. Suppose that An em < 0.1 and that (7-9) hold. Then for m = 
l,...,p-l and j = 0,1,... , sP-"" - 1 

(12) af ) = J2 + ^j?)] [^(1 + ^m)Y 

r=0 

where 

(13) \lm\<ceM + 0{ei^), 
and 

(14) \Vj^U<mdeM + 0{ei^), d = A^-i + (s - 1) c, 
where c is defined by (6). 
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Proof. Let m = 1. Then from (9) it follows that 

s-l 

(15) a« = J2 <^Js+k{l + A^ji) z\ I A« |< As-, 6m, 

A;=0 

which can be rewritten in the following form 

~af = 5] a,-,+,l— [z(l + 71)]'= = J2[ajs+k{l + V^'l) [z{l + 7l)]^ 
where 

From this we obtain 

Now using (11), (15), the definition of d in (14) and the fact that k < s — 1 
we have 

I r/j5 |< (^-1 + {s-l) c) EM + 0{elj) =deM + 0{elj). 

In the same manner, using the equality 

[z{l + 7^)] = [z{l + 7^_i)] (1 + 5m-l)^ 

we get I r^J^ \<mdeM + 0{e^), which is the desired conclusion. ■ 

Theorem 1 Under the assumptions of Lemma 1 the value w{z) computed 
by PEMA satisfies 

N 
n=0 

where 

I /? 1=1 Ip |< ceM + 0{elt), I A„ \< p{As + s c) + 0{elt). 
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Proof. Let m = p — 1. Lemma 1 yields 

^{P-I) ^ ' ^"'[a .,,_,^,(1 + r,%-'^)] [z(l + 7p-l)]^ 

r-=0 

where 

\vt~'^ \<{p-l)deM + 0{e!,). 
By assumptions, we have 

j=0 

This gives immediately the assertion of the theorem. ■ 

So, if the algorithm W satisfies (8), PEMA is numericahy stable in a 
sense (2). 



Table 0: Constants An and Zn for all algorithms. 



Algorithm 


An 


Zn 


Horner 


{c+l)N 





Goertzel 


10Ar2 





PEMA(Horner) 


ps{2c + 1) 


c 


PEMA (Goertzel) 


lOps^ + psc 


c 



N is the degree of the polynomial, c = 1 for real coefficients ao, . . . ,aN and 
c = 1 + \/2 for complex a„. Here p and s are the parameters of PEMA, 
N = s^. Note that for = 2^, the partial polynomials in PEMA are of 
degree 1 and An is of order log2 N, which is a significant improvement 
when compared with the standard versions of both algorithms. 

4 Numerical tests 

This paragraph contains the results of the tests performed in Matlab, ver- 
sion 6.1.0450 (R12.1) with machine precision em ~ 2.2 • 10~^^. We im- 
plemented all methods and compared the results they gave. Of course, it 
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would be the most natural to compare the result given by each of the meth- 
ods with the exact one. However, there are obvious obstructions, i.e. for 
fractional polynomial coefficients or the point z there is no way to obtain 
the exact value of w{z). To deal with these difficulties we used the Matlab 
function "fft", which is perfectly stable (for details see jH], pp. 22-45). The 
function yfft = fft{a) computes Fourier coefficients, namely yk = w{zk), 
k = 0, . . . , N, where w(z) is the polynomial (1) and = ; cj = e is 
the (A^ -|- l)st root of unity: uj^^^ = 1. The values z^ were computed by 
the Direct Call algorithm, i.e. z^ := cos (kt) — ism (kt), t = -^^r^, which is 
known to be very accurate (Cf. [S], pp. 23-24). 
We computed the relative error 

\\y-yfft\\2 

(16 error = 

where y denotes the vector of results given by Horner's, Goertzel's or 
PEMA algorithm for a certain set of points {zj} and yfft is the result 
given by the function "fft" for the same set of points, namely for Zj, where 
j £ {0, 1, 9, 99, 199, 256, 299, 399, 499, 699}. The parameter p in PEMA (see 
section 3) was equal to 2, namely N = s"^ (i.e. s = ^/N). 

The function "fft" can be used provided that l-z] < 1. In general this 
condition is not needed, all algorithms, namely Goertzel's, Horner's and 
both versions of PEMA work for any z £ C 

Figure 1 describes the results for Goertzel's algorithm and PEMA with 
Goertzel's method applied to polynomials with random coefficients. 

Both graphs illustrate the logarithm of error (16) plotted against the 
logarithm of the polynomial degree n, which varies between 2^*^ and 2^^. 
The lower graph represents results given by PEMA, while the upper one 
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Figure 1: Relative errors of Goertzel's and PEMA algorithms 
for polynomials with random coefficients. 



Figure 2 describes similar results for a family of polynomials with coeffi- 
cients given by the formula = f{tk) where = 0.001 A;, k = {),... ,N and 
f{t) = sint + sinlOOt + sinlOOOt. As before the lower and the upper graphs 
represent results given by PEMA and the standard Goertzel's algorithm, 
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respectively. 




Figure 2: Relative errors of Goertzel's and PEMA algorithms 
for polynomials with coefficients Ofc = f(tk) 
where f{t) = sint + sinlOOt + sinWOOt. 



Figure 3 illustrates analogous results for polynomials with coefficients 



Horner and Goertzel 15 
Ojfc = Vk. And again the lower graph represents results given by PEMA. 
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Figure 3: Relative errors of Goertzel's and PEMA algorithms 
for polynomials with coefficients Ofc = Vk. 



Tables 1 — 3 contain values of error (16) for each method and for poly- 
nomials of coefficients given in description above each table. N is the poly- 
nomial degree. The second and the third columns contain results given 
by Horner's rule and the version of PEMA algorithm with Horner's rule, 
respectively. Data in the last two columns is results given by Goertzel's 
algorithm and PEMA with Goertzel's algorithm. This data was used to 
create figures 1 — 3. 

Note that although Goertzel's algorithm gives large errors for large N, 
PEMA using Goertzel's algorithm has much smaller errors; they are com- 
parable with the errors obtained using Horner's algorithm, or PEMA with 
Horner's algorithm. 

Table 1: Relative errors of Goertzel's, Horner's and both versions of PEMA algorithms 
for polynomials with random coefficients. 
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N 


Horner 


PEMA(Horner) 


Goertzel 


PEMA(Goertzel) 


210 


1.6396e - 014 


1.6597e - 014 


6.4827e - 014 


1.6614e - 014 


212 


6.4839e - 015 


6.2312e - 015 


1.7241e - 013 


6.3318e - 015 


214 


6.4597e - 015 


8.8147e - 015 


7.8870e - 013 


8.8450e - 015 


216 


1.0575e - 014 


1.2730e - 014 


1.1884e - Oil 


1.3035e - 014 


218 


3.0060e - 014 


4.3985e - 014 


2.0332e - 010 


4.4917e - 014 


220 


7.1352e - 014 


7.6212e - 014 


5.2591e - 009 


9.7373e - 014 


222 


1.1814e - 013 


1.5060e - 013 


4.1586e - 008 


1.7229e - 013 



Table 2: Relative errors of Goertzel's, Horner's and both versions of PEMA algorithms 
for polynomials with coefficients ak = f{tk) 
where /(t) = sint + sinlQQt + swlOOOt. 





Horner 


PEMA(Horner) 


Goertzel 


PEMA(Goertzel) 


210 


2.1321e - 015 


1.0999e - 014 


4.9016e - 013 


1.1313e-014 


212 


4.3372e - 015 


1.5549e - 014 


3.0108e - 012 


1.7462e - 014 


214 


9.7481e - 015 


2.5365e - 014 


7.0483e - 012 


2.6262e - 014 


216 


3.2760e - 014 


1.0139e - 013 


9.2595e - Oil 


1.1973e - 013 


218 


1.6703e - 014 


1.6408e - 014 


1.0737e - 010 


3.2052e - 014 


220 


8.1798e - 014 


6.0448e - 014 


3.3356e - 009 


8.3002e - 014 


222 


2.6576e - Oil 


3.9179e - Oil 


6.4574e - 006 


4.8041e - Oil 



Table 3: Relative errors of Goertzel's, Horner's and both versions of PEMA algorithms 
for polynomials with coefficients ak = \/k. 
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N 


Horner 


PEMA(Horner) 


Goertzel 


PEMA(Goertzel) 


2io 


5.6281e -015 


5.6566e - 015 


1.5038e - 013 


5.9073e - 015 


2i2 


8.0767e - 015 


8.1583e -015 


3.4601e - 012 


9.3555e - 015 


2i4 


1.8735e - 014 


1.8795e - 014 


1.9228e - Oil 


2.3707e - 014 


2l6 


1.7620e - 013 


4.7930e - 013 


2.1008e -009 


5.2504e - 013 


2l8 


1.1682e -012 


3.5980e - 012 


5.1238e -008 


3.8532e - 012 


220 


8.4972e - 012 


6.1673e -012 


1.4374e -006 


8.1276e -012 


222 


5.1749e - Oil 


4.1890e - Oil 


3.5824e - 005 


5.3874e - Oil 



Appendix. Error analysis of Goertzel's algorithm 

Now we turn our attention to numerical analysis of Goertzel's algorithm. 
Goertzel's method is a special case of Clenshaw's algorithm (Cf. pQ, [2], 
|10)). Our results are similar in spirit to these given by Gentleman who 
gave a floating-point error analysis of Goertzel's algorithm for computing 
Fourier coefficients Yln=o ""^ ^^'^ X^nLi '^n sin with respect to the 

data aQ,...,aj\f and a given argument (Cf. pp. 84-88, [2|). He 

advised to avoid this technique, particularly for low frequencies ^ (e.g. for 
^ = 0). However, we prove that under natural assumptions Goertzel's 
algorithm is numerically stable in a sense (2), as an algebraic polynomial 
evaluation algorithm. These results extend the results obtained in 0, jHl 
for real coefficients a^. Here we consider more general case of complex 
coefficients a„. 

In the exact arithmetic we have for the quantities computed by Goertzel's 
algorithm (Algorithm 2) 

N 

(17) bn = Y.ak\z\'"'' Uk-n{t), n=l,2,...,N, 

k=n 

N N 

(18) u = ^ak\z\'^ Tk{t), V = y ^ak\z\^'^ Uk-i{t), 

fc=0 k=l 
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(19) * G [-1,1] 

I ^ I 

and Tjt(t) and Uk{t) are the Chebyshev polynomials of the first kind and 
the second kind, respectively. They satisfy the recurrence relations (Cf. 

m) 

Tk{t) = 2tTk^i{t) - Tk-2{t), Uk{t) = 2tUk-i{t) - Uk-2{t), k = 2,... 

with To(t) = Uo{t) = 1 and Ti{t) = t, Ui{t) = 2t. 
Moreover, 

Tk{t) =tUk-i{t)-Uk-2{t), k = 2,3,... 
We remind very well known inequalities for t S [—1, 1]: 

(20) \n{t)\<l, \Uk{t)\<k + l, fc = 0,l,... 
It's well known that for |t| < 1 

n{t) = cosk9, c/,(t) = !i^^i^±ll^, 

sm 

where 

t = cose, z = \z\e^\ 6le(0,7r). 

Notice that 

k 

(21) l^=Tk{t)+i^Uk-i{t) for A; = 0,l,...,iV. 

\A \y\ 

Now we analyze numerical behaviour of Goertzel's algorithm in floating- 
point arithmetic. 

Let w{z) = u + iv, bn, p,q denote the quantities computed numerically 
in cfl (see section 3.2). We have 

p = p, q = q{i + i), h\<2eM + o{€lj). 

Therefore, b]\f^i =0, bj\f = a^r and for n = — 1, . . . , 1 we get 

(22) bn = {an + -qn) + Pbn+l + qbn+2, 
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(23) u = {ao + ijo) + xbi + q 62 
where for n = 0, 1, . . . , A/" 

(24) \Vn\<KeMi\an\ + \z\\bn+i\ + \z\^\bn+2\), K = b + 0{eM). 

The constant 5 in (24) is overestimated, but this way error analysis is 
simpler and the essential result is the same. 
Further, we get 

(25) w{z) = {u + iv){l + 52), v = yh{l + 5x), \5x\, \52\ < €m. 
Prom this it follows that 

A' 

(26) In = ^<^k + rik)\zt-'' Uk-n{t), n = l,2,...,iV 

k=n 

and 

N 

(27) u = Y,{ak + r]k)\z\'' Tk{t). 

k=0 

Prom (17), (18) and (20) it follows that 

(28) \u\<9o, \bn\<{N-n + l)gn, n = l,...,N, 
where 

N 

(29) 9n = J2 KIN*^""' n = 0,l,...,7V. 

k=n 

We want to estimate the absolute error \w{z) — w{z)\. Let's write 6„ as 
6„ = 6„ + e„, where from (20), (26) and (27) 

N 

(30) \en\<{N-n + l)Y,\rjk\\z\''. 

k=n 

The formulae (24), (26)- (29) yield 

\r]k\ <KeM (lofel + {N- k)\z\gk+i + {N - k - 1)|^;| W2) + 0{eit). 



Horner and Goertzel 20 
Thus 

(31) \r]k\<2KeMiN-k)gk + 0{eit), for k = 0,l,...,N. 

Now write analogously u = u + eo, where |eo| < X^^o l%l l-^l*^- 

It's easy to check that X^^g 9k l-^l'^ < (-^ + 1) 9o- Prom this and (31) we 

get 

N 

(32) Yl M \z\'' < 2K{N + l)NeMgo + 0{eli). 
Now let's rewrite (25) as 

(33) w{z) = {u + iybi) + ^. 
It's easy to verify that 

(34) \^\<m + l)eM90 + O{el,). 
Further from (21), (26) and (27) we get 

N 

u + iybi = Y {ak + r]k) . 
k=o 

This and (33) yield 

N 



){z)-wiz)\<Y,\m\\z\'' + m. 



1^1 

k=n 

Combining this with (32) and (34) we get the inequality 

\w{z) - w{z)\ < 2K{N + ifeMgo + C(eM), 
which can be reformulated in the following 

Theorem 2 Assume that cfl{z) = z and cfl{an) = an for n = 0, . . . , N . 
Let 

(35) 2K(iV + l)2eM < 0.1, 
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where K is defined in (24). 

Then GoerizeVs algorithm for computing w{z) = Yln=o ^'^ '^^ componen- 
twise backward stable, i.e. 

N 

(36) w{z) = «n(l + ^n) ^ ) l^n 

n=0 

where 

(37) AN = 2K{N + lf.m 

Notice that ^jv ~ lOAT^. Numerical tests in section 4 confirm that the 
constant N"^ is reahstic. 
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